Lecture 8
Doynikov A
February 2, 2016
Анализ Интервенций
Интервенция - резкое внешнее воздействие на динамику временного ряда. Впервые модель интервеций предложена Бокс и Тяо(Box,Tiao) в 1975 году. Посмотрим на график логарифмов авиа миль на пассажира из библиотеки TSA
library(TSA)## Loading required package: leaps## Loading required package: locfit## locfit 1.5-9.1 2013-03-22## Loading required package: mgcv## Loading required package: nlme## This is mgcv 1.8-9. For overview type 'help("mgcv-package")'.## Loading required package: tseries##
## Attaching package: 'TSA'## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tardata(airmiles)
plot(log(airmiles),ylab = 'Log(airmiles)' ,xlab ='Year',lwd=2,col='blue')
Сначала рассмотрим такую модель
\(y_t= m_t+N_t\)
здесь \(m_t\) - функция изменения среднего, \(N_t\) - процесс (взможно ARMA сезонный) без интервенций. Предположим интервеция произшла в момент T. До момента T разумно \(m_t\) положить 0. \(y_t,t <T\) процесс до интервеции. Интервенцию будем строить с помощью функции ступенька
\(S_t(T)=1, t \le T\) и \(S_t(T)=0\) иначе. Функция \(P_t(T)=S_t(T)- S_{t-1}(T)\) называется ипульс. Вместе их называют экзогенными переменными. Если интервенция имела характер мгновенного воздействия, то это моделируется как
\(m_t = \omega S_t(T)\)
, где параметр \(\omega\) подлежит оценке. Если интервенция имеет не мгновенное, а растянутое по времени действие, то это моделируют процессом возможно типа авторегрессии AR(p). Начнем с авторегрессии первого порядка
\[m_t= \delta m_{t-1}+\omega S_{t-1}(T)\] с начальным условием \(m_0=0\). После простейших преобразований это можно переписать как
\[m_t=\omega\frac{1-\delta^{t-T}}{1-\delta}\] при \(t>T\) иначе \(m_t= 0\)
где обычно предполагается, что \(\delta <1\).
В Случае\(\delta =1\) \(m_t=\omega(T-t)\) при \(t>T\) и 0 иначе.
Аналогично с переменной типа ступенька
\(m_t = \omega P_t(T)\)
\(m_t= \delta m_{t-1}+\omega P_{t-1}(T)\) (*)
что после преобразований \(m_t=\omega\delta^{t-T}\) при \(t>=T\)
Удобно использовать оператор сдвига \(Bm_t=m_{t-1}, BP_t(T)=P_{t-1}(T)\)
В частности моель (*) можно переписать так \[m_t(1-\delta B)= \omega BP_{t}(T)\]
или
\[m_t=\frac{\omega B}{(1-\delta B)}P_{t}(T)\]
Также заметим, что \(S_t(T)= \frac{1}{1-B}P_t(T)\) поэтому модель можно записывать как с переменной типа ступень, так и типа импульс.
Перед тем, как приступить к оценке модели интервенции к реальным данным ознакомтесь с возможными моделями в следующем фрейме.
В этом фрейме моделируется модель для \(m_t\) вида
\(m_t=\omega_0I_t(T)+\frac{\omega_1}{1-\omega_2B-\omega_3B^2}I_t(T)\)
где \(I_t(T)=P_t(T)\) или \(I_t(T)=S_t(T)\) В качестве модели для \(N_t\) предлагается использовать AR(1) процесс
\(N_t= \phi N_{t-1}+\epsilon_t\)
Приступим к оценке модели интервенции для ряда airmiles из библиотеки TSA. Предположим, что модель интервенции 11.09.2001 имеет вид
\(m_t= \omega_0 P_t(T)+\frac{\omega_1}{1-\omega_2B}P_t(T)\)
здесь \(T\) это сенябрь 2001 года. Данные имеют сезонную структуру. Выделим данные до терракта и идентифицируем сезонную ARMA модель
wind <- window(log(airmiles),end=c(2001,8))
wind <- log(airmiles)
plot(wind, col = "blue",type = "b",lwd =2)
Хорошо заметен линейный тренд уберем его дифференцированием и построим ACF
dwind <- diff(wind)
plot(dwind, col = "blue",type = "b",lwd =2)
acf(dwind,lag.max=48,lwd =4,col = "blue")
через каждые 12 задерержек (лагов) ACF медленно затухает, что говорит о наличии тренда в сезонности, уберем и его сезонным дифференцированием и опять построим ACF.
sdwind <- diff(dwind,12)
acf(sdwind,lag.max=48,lwd =4,col = "blue")
plot(sdwind, col = "blue",type = "b",lwd =2)
На графике увидим еще интервенции в декабре 1997 года, и декабре 2002, их обсудим позднее. В итоге предлагается оценивать модель \(SARIMA(0,1,1)*(0,1,1)_{12}\) Порядковые номера наблюдений , где заметны интервенции: 12- декабрь 1997,69 - сентябрь 2001,84 - декабрь 2002
air.m1=arimax(log(airmiles),order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12),xtransf=data.frame(I911=1*(seq(airmiles)==69),
I911=1*(seq(airmiles)==69)),transfer=list(c(0,0),c(1,0)),
xreg=data.frame(Dec97=1*(seq(airmiles)==12),
Jan98=1*(seq(airmiles)==13),Dec02=1*(seq(airmiles)==84)),
method='ML')
air.m1##
## Call:
## arimax(x = log(airmiles), order = c(0, 1, 1), seasonal = list(order = c(0, 1,
## 1), period = 12), xreg = data.frame(Dec97 = 1 * (seq(airmiles) == 12), Jan98 = 1 *
## (seq(airmiles) == 13), Dec02 = 1 * (seq(airmiles) == 84)), method = "ML",
## xtransf = data.frame(I911 = 1 * (seq(airmiles) == 69), I911 = 1 * (seq(airmiles) ==
## 69)), transfer = list(c(0, 0), c(1, 0)))
##
## Coefficients:
## ma1 sma1 Dec97 Jan98 Dec02 I911-MA0 I911.1-AR1
## -0.3825 -0.6499 0.0989 -0.0690 0.0810 -0.0949 0.8139
## s.e. 0.0926 0.1189 0.0228 0.0218 0.0202 0.0462 0.0978
## I911.1-MA0
## -0.2715
## s.e. 0.0439
##
## sigma^2 estimated as 0.0006721: log likelihood = 219.99, aic = -423.98Функция arimax методом максимального правдоподобия производит оценку модели вида
\(\Delta^d\Delta_S^D\phi_p(B)\Phi_P(B^s)y_t= m_t+\theta_q(B)\Theta_Q(B^s)\epsilon_t\)
где
\(\phi_p(B)=1-\phi_1B-...-\phi_pB^p\)
\(\Phi_P(B^s)=1-\Phi_1B^s-...-\Phi_PB^{sP}\)
\(\theta_q(B)=1-\theta_1B-...-\theta_qB^p\)
\(\Theta_Q(B^s)=1-\Theta_1B^s-...-\Theta_QB^{sQ}\)
характеристические полиномы авторегресии скользящего среднего, не сезонной и сезонной части соответственно. \(\Delta =1-B\) и \(\Delta_S= 1-B^S\) операторы не сезонной и сезонной разности соответсвено, \(d,D\) - порядки несезонных и сезонных разностей.
В функции arimax параметры \(p,q,d\) и \(P,Q,D\) и период сезонности \(s\) задаются так order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12).
\(m_t\) это матрица столбцы которой \(m_t= (m_{1,t},m_{2,t},m_{3,t},m_{4,t})\) имеют длину равную длины ряда - модели собственно 4 интервенций: сентябрь 2001-порядковый номер в ряде 69;
декабрь 1997-порядковый номер в ряде 12;
январь 1998-порядковый номер в ряде 13;
декабрь 2002-порядковый номер в ряде 84.
Последние три интервенции имеют характер мгновенного вброса, поэтому модель для них выглядит просто
\[m_{k,t}= \beta_kP_{k,t}(T)\] ,\(k=2,3,4; T = 12,13,84\)$
Эти три модели интервенций (вбросов) в функции задаются в R коде параметром xreg=data.frame(Dec97=1(seq(airmiles)==12),Jan98=1(seq(airmiles)==13),Dec02=1*(seq(airmiles)==84))
Матрица \(P_{k,t}(T), k = 2,3,4;\)
xreg=data.frame(Dec97=1*(seq(airmiles)==12),Jan98=1*(seq(airmiles)==13),Dec02=1*(seq(airmiles)==84))
xreg## Dec97 Jan98 Dec02
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 0 0
## 7 0 0 0
## 8 0 0 0
## 9 0 0 0
## 10 0 0 0
## 11 0 0 0
## 12 1 0 0
## 13 0 1 0
## 14 0 0 0
## 15 0 0 0
## 16 0 0 0
## 17 0 0 0
## 18 0 0 0
## 19 0 0 0
## 20 0 0 0
## 21 0 0 0
## 22 0 0 0
## 23 0 0 0
## 24 0 0 0
## 25 0 0 0
## 26 0 0 0
## 27 0 0 0
## 28 0 0 0
## 29 0 0 0
## 30 0 0 0
## 31 0 0 0
## 32 0 0 0
## 33 0 0 0
## 34 0 0 0
## 35 0 0 0
## 36 0 0 0
## 37 0 0 0
## 38 0 0 0
## 39 0 0 0
## 40 0 0 0
## 41 0 0 0
## 42 0 0 0
## 43 0 0 0
## 44 0 0 0
## 45 0 0 0
## 46 0 0 0
## 47 0 0 0
## 48 0 0 0
## 49 0 0 0
## 50 0 0 0
## 51 0 0 0
## 52 0 0 0
## 53 0 0 0
## 54 0 0 0
## 55 0 0 0
## 56 0 0 0
## 57 0 0 0
## 58 0 0 0
## 59 0 0 0
## 60 0 0 0
## 61 0 0 0
## 62 0 0 0
## 63 0 0 0
## 64 0 0 0
## 65 0 0 0
## 66 0 0 0
## 67 0 0 0
## 68 0 0 0
## 69 0 0 0
## 70 0 0 0
## 71 0 0 0
## 72 0 0 0
## 73 0 0 0
## 74 0 0 0
## 75 0 0 0
## 76 0 0 0
## 77 0 0 0
## 78 0 0 0
## 79 0 0 0
## 80 0 0 0
## 81 0 0 0
## 82 0 0 0
## 83 0 0 0
## 84 0 0 1
## 85 0 0 0
## 86 0 0 0
## 87 0 0 0
## 88 0 0 0
## 89 0 0 0
## 90 0 0 0
## 91 0 0 0
## 92 0 0 0
## 93 0 0 0
## 94 0 0 0
## 95 0 0 0
## 96 0 0 0
## 97 0 0 0
## 98 0 0 0
## 99 0 0 0
## 100 0 0 0
## 101 0 0 0
## 102 0 0 0
## 103 0 0 0
## 104 0 0 0
## 105 0 0 0
## 106 0 0 0
## 107 0 0 0
## 108 0 0 0
## 109 0 0 0
## 110 0 0 0
## 111 0 0 0
## 112 0 0 0
## 113 0 0 0Остается наиболее интересная интервенция, сентября 2001 года, ее задают в виде передаточной функции (transfer function) Модель передаточной функции
\(y_t= \frac{u_m(B)}{v_l(B)}x_t+\frac{w_k(B)}{r_h(B)}P_t(T)\)
где \(u_m(B),v_l(B),w_k(B),r_h(B)\) характеристические полиномы степени \(m,l,k,h\) соответственно, \(P_t(T)\) - переменная типа импульс в точке \(T=64\)
\(u_m(B)= u_0-u_1B-...-u_mB^m\)
\(v_l(B)= 1-v_1B-...-v_lB^l\)
\(w_k(B)= w_0-w_1B-...-w_kB^k\)
\(r_h(B)= 1-r_1B-...-r_hB^h\)
Обратите внимание на свободный член в полиномах \(u_m(B),w_k(B)\). Он не равен 1, как мы привыкли видеть в ARMA моделях.
В функции arimax R кода модель передаточной функции задается параметром
xtransf=data.frame(I911=1(seq(airmiles)==69), I911=1(seq(airmiles)==69)),transfer=list(c(0,0),c(1,0)) Так как порядки полиномов заданы \(m =0 ,l=0,k=1,h=0\) То модель передаточной функции выглядеть будет так
\(y_t= u_0x_t+\frac{w_0}{1-r_1B}P_t(T)\)
полагая \(x_t=P_t(T)\),а \(y_t= m_t\)как раз и получим требуемую модель интервенции
\(m_t= \omega_0 P_t(T)+\frac{\omega_1}{1-\omega_2B}P_t(T)\)
Сравним исходные данные с значениями, полученными по модели.
plot(log(airmiles),ylab='Log(airmiles)',col= "blue",lwd =2)
points(fitted(air.m1),col= "red",lwd =2)
Посмотрим на графический фид модели интервенции, которую оценили
Nine11p=1*(seq(airmiles)==69)
inter0911<-Nine11p*(-0.0949)+ filter(Nine11p,filter=.8139,method='recursive', side=1)*
(-0.2715)
plot(inter0911,ylab='9/11 Effects',
type='h',col = "blue",lwd = 2);
abline(h=0)
Ее можно удалить из ряда для разностей
inter0911 = ts(inter0911,frequency = 12,start= c(1997,2))
withoutinter0911= sdwind-inter0911
plot(withoutinter0911,col = "blue",type = "l",lwd = 2)
Данные без интервенции сентября 2001 года показывают, что модель не полностью исключила последствия терракта 9 сентября 2001 года. Но более сложную модель я пока не сделал.